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Given the absence of observations of black hole binaries, it is critical that the full 

q range of accessible parameter space be explored in anticipation of future observation 

with gravitational wave detectors. To this end, we compile the Hamiltonian equa- 
te 

tions of motion describing the conservative dynamics of the most general black hole 
binaries and incorporate an effective treatment of dissipation through gravitational 
radiation, as computed by Will and collaborators. We evolve these equations for 
systems with orbital eccentricity and precessing spins. We find that, while spin- 
spin coupling corrections can destroy constant radius orbits in principle, the effect 



o 



is so small that orbits will reliably tend to quasi-spherical orbits as angular momen- 



turn and energy are lost to gravitational radiation. Still, binaries that are initially 
^ highly eccentric may retain eccentricity as they pass into the detectable bandwidth 

of ground-based gravitational wave detectors. We also show that a useful set of 
natural frequencies for an orbit demonstrating both spin precession and periastron 
precession is comprised of (1) the frequency of angular motion in the orbital plane, 
(2) the frequency of the plane precession, and (3) the frequency of radial oscillations. 
These three natural harmonics shape the observed waveform. 
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I. INTRODUCTION 

Motivated by future gravitational-wave observatories, a campaign to track the inspiral 
of the most generic black hole binaries has long been underway [THTT]. The promise of 
gravitational-wave astronomy lies in our ability to observe and test the full range of astro- 
physical phenomena, including spinning, precessing, eccentric pairs of unequal mass. To 
this end, Will and collaborators have published a series of computations of the equations 
governing black hole binary motion in the Post- Newtonian (PN) expansion to 3.5PN order, 
including spin corrections [8HTT]. We compile those results in an appendix to provide a re- 
source for probing and testing the PN dynamics. We then convert the dissipative terms into 
Hamiltonian variables and suggest a modification of the Hamiltonian equations of motion 
that incorporates the effects of radiation reaction. The modified Hamiltonian formulation 
and the Lagrangian formulation admit equivalent descriptions, though for the purposes of 
investigating the natural harmonics of these systems we will exploit the ease of interpretation 
offered by the Hamiltonian formulation. 

Binary stars that evolve to a pair of black holes can show evidence of eccentricity and spin 
precession in waveforms detectable by LISA [12]. Although long-lived pairs will likely shed 
eccentricity by the time they enter the LIGO bandwidth, the entire orbital plane continues 
to precess along with the spins. Also, black hole pairs formed by tidal capture in globular 
clusters or galactic nuclei may retain significant eccentricity as their signals pass through 
the band of current and future ground-based gravitational- wave observatories [T3J E] . The 
equations of motion of Refs. [8HTT] allow flexibility in handling the gravitational radiation 
emitted by any realistic black hole pair, prior to entering the strong-field. 

There are four questions we can address immediately with this compilation of the inspiral 
equations: (1) Do spinning pairs tend to quasi-spherical orbits? (2) What features are 
generically introduced into waveforms through periastron precession and spin precession? 
(3) How much energy is lost during each burst near periastron passage? (4) How much 
of the orbit and the waveform for eccentric, precessing orbits is well-described by the PN 
expansion? 

The first question (Do spinning pairs tend to quasi-spherical orbits?) must be asked since 
the purely circular orbits are destroyed by spin-spin (SS) couplings. When spin-orbit (SO) 
coupling is incorporated, the entire orbital-plane precesses and the quasi-circular orbits are 
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replaced by quasi-spherical orbits - trajectories that lie on the surface of a sphere whose 
radius shrinks only due to dissipation. However, spin-spin couplings actually destroy even 
these so that there are no constant radius orbits, even if we were to artificially turn off 
dissipation by turning off the half-order terms in the expansion. In other words, if black 
holes spin, there may not be any quasi- spherical orbits and all orbits could retain eccentricity 
at all stages of their inspiral. We can ask how significant the effect is. Since spin-spin is a 
subdominant effect, we find that the eccentric behavior is small and that orbits can appear 
to be very nearly quasi-spherical. 

The second question (What features are generically introduced into waveforms through 
periastron precession and spin precession?) is significant for designing optimal detection 
algorithms and estimating source parameters. For systems possessing spin, modulation due 
to orbital plane precession will leave an imprint on the observed waveform. In the case of 
significant eccentricity, the waveforms are modulated by the radial oscillation of the orbit. 
Characteristically the amplitude is modulated by eccentricity, as is the polarization due 
to the precession of the periastron and the orbital plane. The Fourier transform of the 
waveform will reflect these precessions by reflecting the natural frequencies of the orbit, 
which modulate the frequency evolution from loss of orbital energy that is present in all 
black hole binaries. We demonstrate these features for an example black hole binary that 
possesses all three characteristic frequencies, for which we also address the third question 
posed (How much energy is lost during each burst near periastron passage?). 

The fourth question (How much of the orbit and waveform for eccentric, precessing orbits 
is well described by the PN expansion?) has not been addressed for generic orbits. In the 
limit of quasi-spherical orbits, the dynamics and waveform are usually taken to be accurate 
until the system reaches the innermost stable circular orbit (ISCO), beyond which even the 
conservative dynamics cannot be treated adiabatically. However, the ISCO is formally a 
characteristic of a test-particle orbit in Schwarzschild or Kerr spacetime, and is not well 
defined for binary spacetimes. As we will show, the PN sequence actually diverges well 
outside the ISCO not only for eccentric precessing systems, but for quasi-spherical systems 
as well. We find that the breakdown happens at radial separations r ~ 10M, which is well 
outside the Schwarzschild ISCO (6M), in contrast to the conventional wisdom of using 6M 
as a reference point for truncating PN approximations. 

Due to this limitation of the approximation, the PN equations of motion cannot be used 
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to probe the most extreme form of precession manifest as zoom-whirl behavior - elliptical 
zooms out to apastron followed by multiple nearly circular whirls around periastron [15J. 
Zoom-whirl orbits are most prevalent when periastron drops into the strong-field regime. As 
shown in Ref. [16], complete whirls occur when periastron falls between the IBCO (innermost 
bound circular orbit) and the ISCO. To be clear, zoom-whirl orbits exist and have been 
observed in numerical relativity simulations [17, 18J; they are simply beyond the trusted 
regime for the PN approximation. 

To lay the foundation, we begin with a discussion of conservative black hole dynamics in 
fi7| In §IH[ we discuss the utility of a Hamiltonian formulation of the equations of motion 
for simplifying interpretation of the dynamics of the full inspiral trajectories and resulting 
gravitational waveforms for spinning black hole pairs on eccentric, precessing orbits. In 



§IV[ we explore the limits of our, and indeed any, PN formulation, due to the intrinsic 
non-perturbative nature of the dynamics at the end of the inspiral. In £|Vj we summarize 
the results and the answers to the four questions we posed. For completeness and ease of 
reference, in £jA]the PN corrections to the equations of motion including the dissipation due 
to gravitational radiation are compiled from Refs. [STTTT] as derived by Will and collaborators 
to 3.5PN order with SO and SS coupling. Finally, in ^Bj we convert the PN corrections 
originally computed in the Lagrangian formulation into Hamiltonian coordinates, which is 
the approach we apply throughout the text. 



II. DYNAMICS WITHOUT DISSIPATION 



For this section, radiation reaction has artificially been turned off so the pair shows no 
evidence of dissipation. The advantage of doing so is that we can clearly see intrinsic features 
of the dynamics that dissipation can obscure. In the following section, dissipation is included 
while the insights of this section will continue to guide our perspective. 

If either of the black holes spin, orbital motion famously no longer lies in a plane in the 
absence of special symmetries. The black holes engage in an intricate three-dimensional 
motion tangled by precession. An example of the three-dimensional range of a spinning pair 
is shown in center-of-mass coordinates on the far left of Fig. [TJ We want to highlight the 
two different precession effects: periastron precession and precession of the orbital plane. 
The two effects can be cleanly separated with a natural choice of coordinates [19, .20]. 
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FIG. 1: A black hole pair with milm\ = 1/4. The heavier of the two black holes is spinning at 
0.9 maximal. The lighter is not spinning. Left: The orbit. Middle: The orbit in the orbital plane 
showing the constancy of the periastron and apastron as well as isolating the periastron precession. 
In the orbital plane, the pair precesses around a four- leaf clover. Right: A projection of the orbit 
onto the equatorial plane, showing that the innate symmetry of the orbit is obscured. 

First, consider precession of the orbital plane. If the spins of the black hole and the orbital 
angular momentum are aligned or anti-aligned, motion will lie in the equatorial plane, defined 
as the plane perpendicular to the total angular momentum, J = L + S, where L is the 
orbital angular momentum and S is the sum of the spins. However, when the spins are not 
aligned with the orbital angular momentum, the spins will precess around the total angular 
momentum, and by conservation of the total momentum (i. e. J = 0), the orbital angular 
momentum will precess to compensate (i. e. L = — S). Since the orbital plane is spanned by 
(r, p) and is orthogonal to the orbital angular momentum, L = r x p, as shown in Fig. |2j 
the entire orbital plane precesses. Plane precession is a purely relativistic reflection of black 
hole spins and so isolating the effect draws out distinct signatures for parameter estimation. 

Periastron precession can be isolated, in turn, by following the motion confined to the 
orbital plane. In Ref. [rTHI2U] . it was shown that when there is one effective spin - defined by 
one black hole spinning or two equal-mass black holes with arbitrary spins - that periastron 
and apastron are constants and the periastron precesses at a fixed rate, as shown in the 
middle panel of [Tj By contrast, these precise features are obscured in the equatorial plane, 
as shown by the projection in the far right panel of [T] So, to separate periastron precession 
from the precession of the orbital plane, we select a coordinate system that cleaves the two 
effects. 
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The coordinate system that affects this separation is given by (r, $, where r is the 
radial coordinate, <E> is the angle swept out in the orbital plane, and \I> is the angle swept 
out as L swings around J |IHl ED] • The polar coordinates in the orbital plane that precess 
through space as the plane precesses are (n, $) with 

$ = Lxn . (1) 

where n = r/r. The entire orbital plane then precesses around the total momentum in the 
direction yfr given by 

* i (J x L) 

* = Jx ^ J- . (2) 

|J x L| 

In the set of coordinates (r, $, and their conjugate momenta (P r ,P$,P^), the angular 
conjugate momenta are simply P$ = L and P$ = L z = L • J. The magnitude of the orbital 
angular momentum L is conserved and so, therefore, is P$. 

In the restricted case of one effective spin - again, defined by one black hole spinning 
or two equal mass black holes with arbitrary spins - the component of the orbital angular 
momentum along J, and therefore Pp, is also conserved. The Hamiltonian equations of 
motion then have a remarkably simple form: 

* _ dH p dH 

dP r ' r dr 
dH ■ 

OH 

The natural frequencies given by radial oscillations, f r = 1/T r where T r is the radial period, 
periastron precession, = $/27r, and orbital plane precession fq, = ty/2ir, are functions of 
r only. 

The more general equations, when both black holes spin and have unequal masses, and 
dissipation is included, are collected in the Appendices. In the coordinate system introduced 
above, adding a second spin (for unequal mass) renders = L z no longer conserved, 
although Pj, = L continues to be conserved. Consequently, the precessional frequencies will 
not depend solely on r, but will be modulated by angular position. Adding dissipation drains 
both Pj, and Pj,, although (r, $, remains the natural coordinate system for disentangling 
the two kinds of precession. 



Orbital plane 



FIG. 2: The orbital plane coordinates as defined in Ref. [19J Fig. 1. Left: The orbital plane spanned 
by r x p precesses around the total angular momentum with frequency Right: The coordinates 
as defined within the orbital plane. 

We then convert the radiation-reaction terms derived in the Lagrangian formulation (SjAj 
into Hamiltonian variables (^B|. We are motivated to show our results in the following 
section in Hamiltonian variables by the simplicity of Eqs. (|3| and the simplicity of the 
analytic expressions for the frequencies, as we discuss in the following section (see Eqs. 
([7])). Of course, practitioners are free to choose either the Lagrangian or the Hamiltonian 
formulation, since both admit equivalent descriptions, so that an orbital-plane decomposition 
can be found with corresponding frequencies. The Hamiltonian formulation offers a cleaner 
framework in which to examine the explicit equations, so we favor it slightly for the purposes 
of the paper. 1 



III. GALLERY OF INSPIRALS 

Using the equations of motion of £jBj we can investigate a completely generic orbit. The 
most general scenario allows for both black holes to spin and for those spins to be misaligned. 
These three-dimensional orbits should have three natural frequencies. As described in ^TT], 
the natural frequencies are the radial oscillations of an eccentric orbit 

fr = ^r , (4) 
1 r 

1 For example, the simplicity of L — f x p in the Hamiltonian formulation renders our geometric interpre- 
tation more transparent than the lengthy 3PN corrected definition of L in the Lagrangian variables. 



8 



where T r is the time between successive periastra, the frequency of periastron precession in 
the orbital plane, 

U = 7T , (5) 



and the frequency of plane precession, 



2^ 



(6) 



In the conservative Hamiltonian system, it was found that if spin-spin contributions were 
omitted, then [19] 



2tt/$ = $ = A 



L 



'eff 



L 



2tt/vi> = * 



J X (S e ff X L 

|J x Llr 3 




(7) 



where A = 2dH/dp 2 (Eq. (B17), and S e ff is defined in Eq. (B6). Although we add spin-spin 
coupling, we continue to use Eqs. ([7| as adiabatic estimates for the frequencies. 

The full waveform will be composed of the frequencies Q-(|6]). The waveform is computed 
to leading quadrupole order, although the orbital variables that go into this expression are 
computed in the full 3.5PN system, 



hij 



2n, 
D ' 



ViVj - -TiTj 



(8) 



with D the distance to the sources. This is referred to as the restricted PN waveform 
approximation. For simplicity, we place the Earth along z — J. For greater precision, 
higher-order corrections to the amplitude (including spins) have been computed [ZD, 122] . 
but for our purposes, the leading-order amplitude is sufficient. 



A. Quasi-spherical 

We now consider quasi-spherical, low eccentricity orbits (see also [23])- Low eccentricity 
orbits look quasi-circular in the orbital plane. As the orbital plane precesses, these orbits lie 
on the surface of a sphere whose radius decreases monotonically as gravitational radiation 
is emitted. We work in the Hamiltonian formulation, since the Hamiltonian better lends 
itself to finding constant radius orbits and facilitates the geometric breakdown in terms 
of the orbital plane. Of course, this procedure can readily be repeated in the Lagrangian 
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FIG. 3: The initial conditions are set to those of the constant radius orbit at r, = 50 in the 
absence of spin-spin coupling for an equal mass binary with maximal spins, and spin angles 9\ = 
02 = 45°. Left: The first 10 sec of the orbit is shown (length measured in units of GM/c 2 ~ 
26.5 (M/2OM0)km). Right: The frequency in the orbital plane as defined in Eq. ([7]) marked with 
open circles, and the frequency of plane precession, as defined in Eq. ([7]), marked with crosses. 
Both are shown versus time over the entire inspiral. We plot frequencies in units of (2OM0/M)Hz 
and time in units of (M/2OM0)sec. 



formulation. We find the constant radius orbit according to the prescription in [19] (see also 



). We sketch the argument here. The Hamiltonian of Eqs. (B3) and (B4) does not admit 



a simple effective potential formulation since it is a complicated function of p 2 . Nonetheless, 
we can still use the Hamiltonian as an effective potential at the turning points if we exclude 
spin-spin coupling: 

V cS = H(P r = 0) , (9) 

With this condition, we can find the L of an orbit at a given constant radius. Fixing the initial 
conditions this way and including spin-spin couplings, we can estimate the deviation away 
from quasi-sphericity. We find that the oscillations in the radius due to spin-spin coupling 
are negligible (less than a few percent) and that the radius quickly begins a monotonic 
decrease as gravitational waves are radiated. Therefore, the spin-spin effect seems to be 
too small to induce measurable eccentricity during the inspiral if no eccentricity is present 
initially. 
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FIG. 4: The pair has m^/mi = 1/4, maximal spins, a\ = a2 = 1, and initial spin angles Q\ = 
45°, 62 = 45° and initial values r, = 3000 and = 0.003rj. Left: The first 1150s of the inspiral. 
Right: The final 5s before cutoff. 

As an example, consider an equal mass pair with maximal spins a\ — (1% — 1, and spin 
angles initially misaligned so that 9\ = arccos(L • Si) = 45° and 62 = arccos(L • S2) = 45°. 
This is shown in Fig. |3j where the first ten seconds of the total orbit is shown. The full 
orbit takes ~ 60 seconds to merge for M = 20M Q . We take the initial condition to be 
rj = 50 and set L initially at the quasi-spherical value determined by Eq. (|9]). The natural 
frequencies /$ and are shown in the right panel. Since the orbit is quasi-spherical, the 
frequency of radial oscillations f r is not informative. Periastron precession dominates over 
plane precession, which is expected since plane precession is a higher-order effect. 

Using physical units, with the initial separation of ~ 50, the gravitational radiation 
has an initial frequency 

2x/ $ ~ 9.5Hz (10) 

which is already nearing the Advanced LIGO/VIRGO band for two 1OM black holes, and 
is sweeping through the bandwidth as the orbital separation decreases. The much slower 
frequency of precession, does not have a significant effect on the waveform for the quasi- 
spherical case. 
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FIG. 5: Strain for an optimally-oriented binary with M = 20Mq at D = 100 Mpc. From left to 
right, we show three periastron passes occurring ~ 3000 sec prior to waveform truncation (shifted 
in time by ~ 34 sec each to increase their overlap), the final 10 seconds, and the final ~ 3000 
seconds of the waveform. 
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FIG. 6: AH is the absolute value of the difference between the energy as a function of time and 
the initial energy. Left: The burst of energy emitted during one periastron passage (r p ~ 38). As 
the pair separate to apastron r a ~ 3000, negligible energy is lost as confirmed by the flatness of 
AH. Each burst is small although a few percent of the total energy in [x is lost cumulatively prior 
to cutoff as shown on the right. 



B. Eccentric, unequal masses 

We next consider an eccentric pair with mass ratio mi/1712 = 1/4. Each black hole is 
spinning maximally. Both spins are misaligned with the orbital plane with 9\ = arccos(L • 
Si) = 45° and 8 2 = arccos(L ■ S2) = 45°. The pair begins at large separation, Tj = 3000, 
and with angular momentum L = 0.003?v The first periastron passage is r p ~ 38, so the 
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FIG. 7: Left: The three natural frequencies versus time for the eccentric configuration. There are 
two lines for one marking the value at periastron (top line) and the other marking the value 
at apastron (bottom line). Same for fa. All frequencies are plotted in units of (20M Q /M)Hz and 
all times in units of (Af/2OM0)sec. Right: The characteristic strain, h c , of an optimally-oriented 
realization of the eccentric configuration, at a distance of 100 Mpc, again scaled for M = 20Mq. 
We show h c for each of the first three periastron passages, the final 100 sec, and the entire 3000 
sec intervals shown in Fig. [5j We show the LIGO and Advanced LIGO sensitivity curves |24j for 
comparison. 

pair is set on a highly eccentric orbit. The eccentricity can be denned loosely as 



where r a is the apastron and r p is the periastron. The eccentricity begins large, e = 0.97, 
and drifts down to e = 0.17 by the time the simulation is ended. Cutoff is set when the 
3.5PN terms become larger than the 2.5PN terms, or, equivalently, when E > 0, as we 
discuss in the next section. At cutoff, the periastron is r p ~ 12 and the apastron is r a ~ 17. 
Our choice of the initial pericenter distance and eccentricity is consistent with the capture 
formation scenario of stellar mass black hole binaries in globular clusters [14] . Therefore, 
an initially highly eccentric orbit could retain significant eccentricity by the time it evolves 
into the bandwidth of a terrestrial network of gravitational-wave interferometers. 

The pair experiences both periastron precession and orbital-plane precession during a 
total of 1384 orbits, some of which are shown in Fig. [4j On the left of Fig. [4] the first 
1150s of the orbit are shown and on the right the final 5s of the orbit is shown before the 
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simulation is cutoff. Both the precession of the periastron and the precession of the orbital 
plane are apparent in Fig. |4| The precession of the orbital plane is more prominent the more 
disparate the masses of the black holes. While spin precession remains substantial to the 
end, periastron precession becomes less significant as eccentricity is shed. 

Since the orbit is highly eccentric, a burst of gravitational radiation is emitted near 
periastron passage as shown in Fig. |5j During the most eccentric cycles, the bursts are what 
one might expect from a parabolic encounter, with energy and angular momentum losses 
occurring primarily near periastron. A tiny amount of energy is emitted per periastron 
passage for the M = 2OM pair shown in Fig. |6j although a few percent of the total /i is 
shed prior to plunge. The bursty waveform is highly distinct from the extremely regular 
quasi-circular waveforms, clearly showing effects of the precession. 

The natural frequencies /$ > f r > are shown on the left of Fig. [7j Both and fa 
are shown only at periastron (top line) and apastron (bottom). Therefore, there are two 
lines each for and fa. The value of the frequencies at apastron is much smaller than 
at periastron if the eccentricity is large. These discrete values approach each other as the 
eccentricity decreases. Since T r is defined as the time between successive perihelia, there is 
only one line for f r = 1/T r . The wobbling in fq, is a result of spin-spin coupling. The time 
between bursts is given by the radial period, T r = l/f r , and the frequency of the burst is 
given by fq> at periastron passage during the highly-eccentric, nearly parabolic encounter, 
and by 2 x when the orbit is less eccentric. 

The right panel of Fig. [7] shows the characteristic strain, 



h c = 2f 



2 



h. 



2 



(12) 



for the three time intervals shown in Fig. |5j as well as the noise amplitude spectral densities 
for the initial and Advanced LIGO detectors. The self-similarity of the periastron passages 
is evident in Fig. [7} We suggest the possibility that this self-similarity could be employed to 
coherently combine the signals from all the periastron passages of a highly eccentric system, 
while ignoring the data in between where the signal power is negligible. 

The characteristic strain of the full signal has a flat plateau that spans the frequency 
range between the spin precession frequency and twice the orbital frequency, both evaluated 
at periastron, because the power emitted during each orbit is dominated by the emission 
at periastron. The characteristic strain of the final 10 sec shows a very clear step structure 
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due to the large set of harmonics excited by the eccentric motion. We emphasize that 
these harmonics are not the instantaneous harmonics from PN corrections, as we have not 
included higher-harmonic PN corrections. Rather, the harmonics apparent in Fig. [7] are the 
result of the large change in instantaneous frequency over the course of each orbit, which 
is represented in Fourier space as higher (and lower) harmonics of twice the mean orbital 
frequency. Again, the highest plateau spans between the spin precession frequency and twice 
the orbital frequency, with the step at lower frequency occurring at the orbital frequency, 
and the many steps at higher frequencies occurring at larger integer multiples of the orbital 
frequency. 

IV. DOMAIN OF VALIDITY FOR THE PN APPROXIMATION 

The PN equations of motion are a series expansion of corrections to the acceleration 
of the binary components, with whole order terms representing conservative, relativistic 
corrections, and half order terms representing dissipative corrections. As such, the expansion 
should not be trusted when higher-order terms become comparable in magnitude to lower- 
order terms. In this regard, it is noteworthy that the sign of the acceleration term along 
the radial direction is negative at 2.5PN order, as expected, and the black holes are drawn 
together. However, the sign of the acceleration term along the radial direction is positive 
at 3.5PN order. In other words, the 2.5PN terms over estimate the effects of dissipation, 
and the 3.5PN terms temper this over estimate. A consequence of the breakdown of the 
PN expansion in the strong-field regime is that the 3.5PN term comes to dominate, and 
spuriously drives the system to a larger separation, with the 3.5PN term acting as a source, 
rather than a sink, of energy and angular momentum. Implementation of the equations of 
motion must be cut off before this happens, as the PN approximation has already broken 
down at this point. We find that the break down tends to happen at radial separations 
r ~ 10 in units of the total mass, which is well outside the Schwarzschild ISCO at 6M that 
is frequently used to represent the boundary for a valid application of the PN expansion. 

All of the simulations in the preceding section are therefore cut off when the PN approx- 
imation begins to break down. Specifically, we use the criterion that the simulation ends 
when the rate of energy loss drops below some threshold, indicating that the 3.5PN correc- 
tion is approaching the same magnitude as the 2.5PN term. In Fig. |HJ we show the relative 
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magnitude of different PN-order contributions to the (Lagrangian) equations of motion. We 
focus on two cases: an equal mass, nonspinning binary undergoing quasi-circular inspiral, 
and our eccentric configuration from the previous section. The first case is likely to be a 
best-case scenario for the convergence of the PN expansion, while the latter is likely to over 
extend the PN expansion at larger orbital separations. We find that, for the conservative 
terms, the ratio of adjacent whole-order terms in the PN sequence remains less than unity 
over the full domain, and indeed would remain so down to the aforementioned ISCO radius. 
However, the ratio of the two half order, dissipative terms exceeds unity at r = 9.125M in 
the circular case, and at r = 10M in the eccentric spinning case. Since the approximation is 
certainly unreliable prior to the point where this ratio equals unity, we suggest that generally 
the PN expansion cannot be reliably applied for separations r < 10M for any configuration 
when dissipation is included. 

Another possibility is that the optimal asymptotic expansion for the PN sequence should 
be truncated at 3 PN order. This is an interesting possibility, and would in fact be consistent 
with findings in the Schwarzschild test-mass limit, where much higher PN terms are available 
to demonstrate this behavior |25j. Since for comparable masses we lack higher PN-order 
terms, it is not yet possible to distinguish between these two possibilities. 

V. SUMMARY 

We have compiled the equations of motion of Refs. [8HTT]. and have transformed them into 
Hamiltonian coordinates for ready-to- use equations governing general black hole binaries. 
With these in hand, we can study any black hole pair, including the effects of spin, eccen- 
tricity, and their accompanying precessions. Given this flexible system, we have answered 
four basic questions: 

• 'Do spinning pairs tend to quasi-spherical orbits?' 

We find that the deviation for a constant radius orbit due to spin-spin coupling is 
insignificantly small, so quasi-spherical orbits are accessible. Consequently spinning 
pairs do tend to quasi-spherical orbits. 

• 'What features are generically introduced into waveforms through periastron and spin 
precession?' 



16 



1 PN/O PN, circular 




r [M] 

FIG. 8: The ratio of different PN-order contributions to the Lagrangian equations of motion for 
a circular, equal mass, nonspinning inspiral, and separately for our eccentric configuration. Even 
for the circular case, the 3.5 PN dissipation term grows larger than the 2.5 PN dissipation term 
for radial separations less than ~ 9M. Less ideal cases, not surprisingly, appear to make the PN 
sequence diverge at larger radii. 

As was shown in [19], an orbit can be projected onto an orbital plane. The char- 
acteristic frequency within the plane supplants the usual coordinate frequency. The 
frequency of plane precession and the frequency of radial oscillations provide the other 
crucial frequencies for analysis. We have shown that these three frequencies leave a 
characteristic imprint on the waveforms in Fourier space. Furthermore, the frequency 
of bursts for highly eccentric orbits and the quiescent time between bursts directly 
reflect these natural harmonics and could encourage novel data analysis techniques. 

• 'How much energy is lost during each burst near periastron passage?' 

Given the caveat that most of the energy will be lost during the merger, which is 
beyond the reach of the PN approximation, for highly eccentric binaries we show that 
less than 10~ 3 of a percent of the reduced mass /i is lost to gravitational radiation per 
periastron passage. Some pairs may execute thousands of orbits before plunge so that 
a few percent of /i is lost during inspiral up to r ~ 10M. 

• 'How much of the orbit and waveform for eccentric, precessing orbits is well described 
by the PN expansion?' 
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The PN expansion has broken down whenever higher-order corrections become larger 
than lower-order corrections. We find this tends to happen at r ~ 10M, which is well 
outside of the Schwarzschild ISCO that is often used to demarcate the breaking point 
of the PN approximation. Possibly the PN sequence diverges inside ~ 10M, or it may 
be that the optimal PN expansion inside this radius should be truncated at no higher 
than 3 PN order. 

The orbits shown give a sample of the complete range that can be probed in detail for 
generic, spinning, precessing, eccentric black hole pairs. We intend to build a data set 
of black hole pairs to be made available as a testing ground for developing data analysis 
techniques. 
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Appendix A: Lagrangian Equations of Motion 

In this section we compile the equations of motion including spin couplings and dissipation 
as computed by Will and collaborators over the series of papers [8HTT]. We measure length 
in units of total mass M = mi + and write the equations in terms of the dimensionless 
center-of-mass coordinate x and center-of-mass velocity v. The reduced mass is defined as 
7] = /j,/M = m\m2jM 2 . The Lagrangian equations of motion in x, v with spins added are 
(with x = a/x ■ x the harmonic radial coordinate and x = n ■ v) as compiled from Refs. 
El EH [22] 



x = 

X = 



(Al) 

a N + &ipn + &2PN + a 2.5P7V + a 3PV + a 3.5P7V + &PN-SO + a 3.5PN-SO + &PN-SS + a 3.5PN-SS 



and 



n 

a^v = ~ 

ar 



aipjv = 1(1 + 3^)v 2 — 2(2 + 77)- — Ivx 2 } + ^2(2 - V )x 
x 2 I x 2 ) x 2 

*2PN = ~ 1 1 (12 + 29^)^ + 77(3 - 477) (v 2 ) 2 + " St?)* 4 

3 1 v 2 A 2 
- -77(3 - 4r/)v 2 i 2 - -77(13 -477) (2 + 25r/ + 2r/ 2 ) — 



via; 
+ ^ 1 2 



77(15 + 4r/)v 2 - (4 + 41r/ + 8r/ 2 )- - 3r/(3 + 2r])x 2 

x 



n f8 /1\ /171 

1 5^ V a; / V a; 



3- + w 2 



— 



+ 



n 



ar 



16 + 



1399 41 



12 16 



71 2 
vr" I 77 + —77 



20827 123 

1 7T 

840 64 



1 + 



/22717 615 



I" 



168 



64 



11 



H -it 77 H 77^ - 77? d 



x 



X 



+ | (11 - 4977 + 52T7 2 ) v 6 



35 



, 77(1 - 577 + 5t7 2 )x 6 + I (75 + 3277 - 4O77 2 ) 



+ | (158 - 6977 - 6O77 2 ) Q i 4 - 77 (121 - I677 - 2O77 2 ) Q 

-^77 (20 - 7977 + 6O77 2 ) 7; 4 i 2 + ^77 (4 — I877 + I777 2 ) v 2 x 4 1 
8 '8 J 

G) 2+ i( 65 - 152 "- 



V 

H — ~x 



x 



15 



5849 123 



840 32 



tt z 77 - 25t7 z - 8t/ j 



77 (3 - 877 - 2t7 2 ) x a + 77 (15 + 2777 + IO77 2 ) Q t; 2 
(329 + 17777 + IO877 2 ) Q A 2 - ^77 (16 - 3777 - I677 2 ) v 2 x 
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&3.5PN 



fx I 8 (I 

: T] I — 1 x 



x^ 



X 



23 / 1 \ 3 

— (43 + Urj) l-j + — (61 + 70t/)t; 4 + 70x 4 



+ — (519 - 1267?7) ( - )v 2 + - (147 + 18817) ( - ] i 2 - — (19 + 2r/) v 2 x 2 
42 V x J 4 V x / 4 



+ 



I / 1 \ I 

— (1325 + 546r7) l-j + — (313 + 42r/) v 4 + 75x 4 



-^(205 + 777r/) Q |tr 



— (205 + 424r/) ( - ] x 2 - - (113 + 2r/) d 2 x 2 
12 V x / 4 



(A2) 



using the 3PN and half-order terms from Ref. [9]. 

We also absorb an M 2 into the spins so the physical spin in the Lagrangian system is 
denoted Sj = ai{m 2 /M 2 ) where Oi is the dimensionless amplitude < |aj| < 1. As we will 
see in £jBj this is a different normalization from that for spins in the Hamiltonian case, where 
spins are expressed as Sj = ai(m 2 j\iM). The spin-orbit contribution to the acceleration is 

ap JV -so = ^|6-L J v(S + 0-vx (4S + 3£) + 3xn x (2S + 

X I X 

with £ = (m2 /mi)§i + (mi/ 7712)^2 [10J. It is customary to define a reduced Newtonian orbital 
angular momentum, 

Ljy = x x v . (A3) 

The spin-orbit dissipation term with the covariant spin supplementary condition is (see the 
appendix of Ref. [TU] ) , 



&3.5PN-SO 



T] I xn 



P» <-y»4 I rf 

.2 



^120t; 2 + 280x 2 + 453-^ ) L N ■ S 

r 



+ ( 120?/ + 280± a + 458- ) L x • < 



+ - I 87v 2 - 675x 2 - — - ] L^v • S + 4 ( 18v 2 - 150x 2 - 66- J L^v ■ f 
1 [\ 3 x J \ x J 

-^iv x S ( AZv 2 + 15x 2 + 364- J + Jivx^ ( 291i; 2 - 705x 2 - 772- J 



+ x § ( 31c 1 - 2(>0r 2 ;- 2 -- 2 lar 1 



689 ,1 



v A - + 537x z - + 



1 4 1 



+ x £ ^115tj 4 - 1130v 2 x 2 + 1295x 4 - ^V- + 849x 2 - + 



3 x 
869 J 
3 x 



1 ' 3x 2 
1 44 1 

x ' 3 x 2 



(A4) 
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And the spin-spin contribution is 
3 



&PN-SS 



n(Si • S 2 ) + Si(n • S 2 ) + S 2 (n • Si) - 5n(n • Si)(n • S 2 )] (A5) 



See also [22] • The 3.5PN order effects of Spin-Spin coupling from Ref. [H] are 



a 3.5PN-SS 



X° 



11 



287x 2 - 99^ 



( 



1029A 2 - 123v 2 



541 1 
5 x 
629 1 
lOx 



171 



195x 2 - 67- (Si • S 2 



x 



J x(Si ■ S 2 ) - ^2646A 2 - 7Uv 2 + i(n • S 1 )(fi ■ S 2 ) 

((n • Si)(v • S 2 ) + (n • S 2 )(v • Si)) - 336i(v • Si)(v • S 2 ) 

1038 1' 



174t> 2 - 1386x 2 



5 x 



[h ■ Si)(n • S 2 ) 



-438A ((n ■ Si)(v ■ S 2 ) + (n ■ S 2 )(v ■ Si)) + 96(v • Si)(v • S 2 )] 



+ 



27 
— ? 
10 



75 . 2 509 1 
~2 X ~ l30"x 

-x 2 



((v • S 2 )Si + (v • Si)S 2 ) 
x ((n • S 2 )Si + (n • Si)S 2 ) 



(A6) 



2 2 10 x / 

The spin supplementary condition has no effect on spin-spin terms up to this order. We have 
not yet included the quadrupole-monopole contribution. Finally, the spins precess according 
to 



Si 
S 2 



(x x v) x Si / 3m 2 
3mi 



x 



(x x v) x S 2 



2 + 



2m 2 



We could also add to the right hand side of (A7) the spin-spin terms: 

1 



(§i)pn-ss = * (S 2 - 3(n • S 2 )n) x S x 



x-' 



1221 and 



(S 



l)3.5PN-SS 



X° 



(v • S 2 ) + 30x(n • S 2 ) (n x Si 



(A7) 



(A8) 



(A9) 



[TT] . Equations (A1)-(A9) constitute the Lagrangian dynamical system. These equations 
are complete through 3.5PN order except for the spin contributions, for which additional 
terms have been calculated (see for instance Ref. [2S1 12Z])- Since spin effects are already 
small at the order we include here, we have not pursued inclusion of the higher- order spin 
terms. 

Notice these spins are not reduced by [i and a definition of dimensionful angular momen- 
tum will have a /x in it: J = fiL^ + ... + Si + S 2 . 

Next we re-express the dissipation terms in the language of the Hamiltonian system. 
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Appendix B: Hamiltonian Equations of Motion: Including Dissipation 

The Hamiltonian PN-formulation is by definition conservative and therefore does not 
incorporate dissipation [THO]. We can, however, take the half-order acceleration terms from 
the Lagrangian coordinates and simply convert them to the coordinates appropriate for the 
Hamiltonian. We begin by first laying out the whole-order terms in the usual Hamiltonian 
system p~H6]. 

In a Hamiltonian formulation, the equations of motion are derived from 

r = — p = -— (Bl) 
dp ' dr 

As is standard convention, we work in dimensionless coordinates: the dimensionless coor- 
dinate vector, r, is measured in units of total mass, M = mi + m 2 , for a pair with black 
hole masses m% and m 2 . The canonical momentum, p, is measured in units of the reduced 
mass, \x = miiri2/M. The dimensionless combination rj = fi/M will again prove useful. We 
write vector quantities in bold. The coordinate r is to be understood as the magnitude 
r = y/r ■ r. Unit vectors such as n = r/r will additionally carry a hat. Finally, we have 



used the dimensionless reduced Hamiltonian if = in Eqs. (Bl ), where Ti is the physical 



Hamiltonian, to 3PN order plus spin-orbit terms [THB]. H can be expanded as 

H = H N + H 1PN + H 2 pn + H 3PN + H so + H S s , (B2) 



where 
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H N = 



(B3) 



Hipn = I (377 - 1) (p 2 ) 2 - ± [(3 + 77) p 2 + r)(n • p) 2 ] J + ^ 
^2Piv = ^ (1 - St? + St? 2 ) (p 2 ) 3 + 11(5- 20 V - St? 2 ) (p 2 ) 2 



-2rf(n • p)V - 377 2 (n • p) 4 ] ± 

+ + 877) p 2 + 377(n • p) 2 ] 1 - J (1 + 377) 1 



#3PW = 7^7 (-5 + 35r/ - 7O77 2 + 357/ 3 ) (p 2 ) 4 + — (-7 + 42t? - 53t7 2 - 5t7 3 ) (p 2 ) 3 



1 

128 



+ (2 - 377)77 2 (n ■ p) 2 (p 2 ) 2 + 3(1 - ^ 2 (n ■ p)V - 577 3 (n ■ p) 6 ] 
+ 



l(-27 + 13677 + 10977 2 )(p 2 ) 2 + 1(17 + 3077)77(11 • p) 2 p 2 + 1(5 + 4377)77(11 • p) 4 



+ |^1 h 600 + ( 37r ' " 1340 ) V ~ 552T7 2 ] P 2 - 1 (340 + 3tt 2 + II277) 77(11 • p) 2 | 1 



H 



so 



+ 1 [12 + (872 - 63vr 2 ) 77] 1 

L • S e ff 



(B4) 



and adding spin-spin, 



H. 



ss 



[ 3(Si-n)(S 2 -fi) -Si-S 2 + 



m 2 



mi 



(3(Si • n)(Si • n) - Si • Si) + — ^- (3(S 2 • n)(S 2 • n) - S 2 • S 2 ) 

Notice that H ss includes Sf and Sf terms. For two spinning black holes S e ff is 2 



(B5) 



S e ff — ^lSi + ^S 2 



(B6) 



where the dimensionless reduced spins are defined as 



Si = ai(m 2 //iM) , S 2 = a 2 (ml/fiM) 



(B7) 



and 



. ( 3m 2 



3777.1 



(B8) 



2 The definitions for S e fr can vary in the literature up to an overall constant although the reduced Hso 
must be the same for all prescriptions. 
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The dimensionless spin amplitudes are confined to the range < ai )2 < 1. The spins precess 
according to 



rigi _ dH 



x Si 



(B9) 



or 



S 1 = 5 1 



L x Si ji 



■ L x S 2 n 

&2 = 02 5 1 o 



3n ! ( S 2 + — Si ) • n 

mi 



3n( (Si + ^S 2 ) n] - S, 

m 2 



The orbital angular momentum precesses according to 



x Si 
x S 2 



L=rxp+rxp 

L X S e ff \i 



3fi 



S 2 H Si I • n 

mi 



x Si 



3n 



Si + ^S 2 

m 2 



n 



(BIO) 

(Bll) 
x S 2 



Adding these together, it follows that the total angular momentum J = L + Si + S 2 is 
conserved in the absence of dissipation so that the orbital angular momentum precesses 
according to 

L = -Si - S a • (B12) 

Now, we are ready to convert the radiation-reaction terms of appendix [A] into Hamiltonian 
variables. To do so, we need to relate the Hamiltonian variables (r, p) to the Lagrangian 
variables (x, v). To convert the 2.5PN radiation-reaction term, we only need the lPN-order 
coordinate conversion to catch all corrections up to 3.5PN. To convert the 3.5PN radiation- 
reaction term, we only need the zeroth-order PN coordinate conversion. It has been shown 
that to 1PN order [? ] 

r = x , (B13) 

so that 

r = x = v . (B14) 

What we really want are the harmonic variables of the Lagrangian formulation (x, v) in 
terms of canonical Hamiltonian variables (r, p) to 1PN order. We have x(r, p) above. To 
find v(r,p), we use Hamilton's equations: 

S eff x r 



r = Ap + Bh 
p = Cp + Dh 



S eff x p 3n dH ss 

O 1 tlSO 7, 

r 6 r or 



(B15) 
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where A,B,C,D are 



.4 
B 
C 
D 



n dH PN 




dp 2 


r,(n 


dH PN 




d(h ■ p) 


r.p 


1 dH PN 


r d(h ■ 


P) 


dH PN 




dr 


P,( 


dH PN 




dr 


P.( 



(B16) 



r.p 



+ 



B 

r 

dH PN 



d(h ■ p) 



(np 



r.p 



(fi-p)C • 



(B17) 



and Hpn = Hn + Hipn + H 2 pn + H^pn- Explicit expressions can be found in [20J, but we 
will only need A and B here to 1PN order to calculate v(r, p) to 1PN order and ultimately 
find the radiative contributions to the accelerations. Taking the appropriate derivatives of 
the Hamiltonian to 1PN order, we have 



A <1PN = 1 + \ (3r7 - 1) p 2 - (3 + 77) - 
2 r 



B 



<1PN 



-77(n-pj 



This gives us v in terms of (r, p): 



1 1 \ 1 

1 + o ( 3r l ~ X ) P 2 ~~ ( 3 + V) ~ P - V (n • p) -n + 

2 r J r 



Seff x r 



(B18) 



(B19) 



and 



1 + (3?7 - 1) p 2 - 2 (3 + 77) 



P 



2 n /« x2 1 , P • (Seff X r) 

277 (n-p) - + 2- 



r 



1 + (3r/ - 1) p 2 - 2 (3 + 77) - p 2 - 2r/ (n • p) 2 - + 2 



'eff 



We can now re-write a 2 5 p7v of Eq. (A2) in terms of ADM variables as 



a 2.5PN^ADM 



5 Vr, 



(B20) 



(B21) 



11 



17 n , \ a , , p (n • p) S e ff • L 

- + 3p 2 + 3 377 - 1 p 4 - 6 3 + 77 — - 6r/^ ^- + 



P 



- + P 2 + |(37/-1)P 4 + J(377 - 21)^ 
r 2 2 r 



3(3 + ,) '-a,^ 



77(n-pJ 

2 Seff • L 



25 



We can also replace r with 



r = n ■ v 



1 1\ 1 

1 + g ( 3r ? ~ X ) P 2 ~~ ( 3 + ^ ~ J fi ' P ~ ^ ( ft ' P) ~ 

1 + ^ (3^7 - 1) P 2 - (3 + 277) - )n-p 



(B22) 



To get 



&2.5PN^ADM 



5 Vr 



n 



P 



np 

3 

- + 
r 



17 

3r 



9 



+ 3p 2 + - (3?7 - 1) p 4 - - \hr) + 



179 \ p 2 1 



(B23) 

^2 



- - (25?7 + 51) 



6rj 



v n 'PJ 



P 2 + o - X ) P 4 + o ( 3r ? - 21 



p^ 



3(3 + 7;)- -277 



n- PJ 



>eff 



S e ff x r / 3 



Since we used a coordinate change valid to 1PN order, Eqs. (B23 and B21) includes 3.5PN 
corrections as well as 2.5PN corrections. 

To rewrite sl3.5pn^adm we simply take the harmonic expressions and replace x = r and 
v = p, noting that L^r = L to lowest order, and that the spins of Refs. pHl E] are rj times 
those of this section: Sj = T]Si. 

Now, to incorporate the effects of dissipation, we can add the appropriately converted 



accelerations. Taking the time derivative of Eq. (B19) and rearranging to solve for p, there 
are additional corrections of the form 

p = v ( 1 - \ (377 - 1) p 2 + (3 + rj) -) - (377 - l)(p • p)p + 77(11 • p)-n + ... (B24) 
\ 2 r J r 

We then use v = a2.5p7v_5.ADA/ + &3.5Pn^adm and drop all terms higher than 3.5PN to get 

p = ... + a. 2 . 5 PN->ADM ( 1 - g ( 3r ? - 1) p 2 + (3 + 77) - J + &3.5PN-+ADM+ 



(377 - l)(p • a 2 . 5P7 v)p + 7/(n • a 2 . 5 pjv)-n 

r 



(B25) 



where we have only isolated the relevant half-order terms so ..." represents all the whole- 



order PN terms. Here ^.^pn^adm is defined in Eq. (B23), and we replace p where it 



occurs on the right hand side of (B25) with a 25 p/v evaluated at x = r, v = p. The term 



&3.5PN-^adm is also evaluated at x = r, v = p, and, again, the spins of the previous section 
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are related to these spins through Sj = rjSi, yielding 



a : v,/> \ ... \v \i aL 2 .5PN->ADM ( 1 - ^ (3i7 - 1) p 2 + (3 + 17) ^ ) = 
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5 Vr 



x 













np | 


[r 2 





- + 3p +3(377-l)p -2(4r ] + 9) v --^-Qr ] —— + Q— r 
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- + p 2 + (377 - 1) p 4 - 2 (77 + 3) — - 277 

r r r 



+ 2 



r 3 r 2 

S e ff ' L 



S eff x r / 3 



Finally, our equations of motion become 

A . n - . S eff x r 
r = Ap + Bn H — — 

. ~ S e ff x p 3n dH ss _ 

p — Op + iJn H 1 U S 7\ r &-2.5PN—tADM + &3.5PN->ADM 

r 6 r or 

+ 77(11 ■ a 2 . 5 p7v)^n - (377 - l)(p • a 2 . 5 pjv)p 

We can regroup the accelerations to 2.5 order and 3.5 orders in ADM variables: 

A . D . . S efr x r 
r = Ap + Bn H — — 

n . n ~ . S eff x p 3n dHss . (2.5) . (3.5) 

p = Cp + Dn + — - — + — H so ~ -?^r + a-ADM + & adm 



(B26) 



(B27) 



with 



and using 



(2.5) _&(V\ fn 
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n p 



17 
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+ 3p< 



(B28) 
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877 n p 2 

3r^7 
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(B30) 



and the 3.5PN piece in 5 2 .5Pat->.adm 
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1 r 2 
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(B31) 
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and adding all of these 3.5PN pieces to 
1 



a3.5PAT-s.ADM — rV 
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Then there are the spin pieces: 
.(3.5PJV-50) _ _ V 2 /(n-p)n [/ 
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and 
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(B33) 



with ^adm = (rn, 2 /mi)Si + (mi/m 2 )S 2 . 

Given the equations of this section and the previous section, the equations of motion can 
be evolved in either set of variables. 



In ^ITj we introduced the (r, $, \l/) coordinates. The equations of motion (B28) can be 
projected onto this basis as was done in Refs. [191 120] - 
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